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The random phase approximation to the correlation energy often yields highly accurate results 
for condensed matter systems. However, ways how to improve its accuracy are being sought and 
here we explore the relevance of singles contributions for prototypical solid state systems. We 
set out with a derivation of the random phase approximation using the adiabatic connection and 
fluctuation dissipation theorem, but contrary to the most commonly used derivation, the density is 
allowed to vary along the coupling constant integral. This yields results closely paralleling standard 
perturbation theory. We re-derive the standard singles of Gorling-Levy perturbation theory [Gorling 
and Levy, Phys. Rev. A 50 , 196 (1994)], highlight the analogy of our expression to the renormalized 
singles introduced by Ren and coworkers [Ren, Tkatchenko, Rinke, and SchefHer, Phys. Rev. Lett. 
106 , 153003 (2011)], and introduce a new approximation for the singles using the density matrix 
in the random phase approximation. We discuss the physical relevance and importance of singles 
alongside illustrative examples of simple weakly bonded systems, including rare gas solids (Ne, Ar, 
Xe), ice, adsorption of water on NaCl, and solid benzene. The effect of singles on covalently and 
metallically bonded systems is also discussed. 

PACS numbers: 71.15.-m, 71.15.Nc., 71.15.Dx, 71.55 Gs 


I. INTRODUCTION 

In the last decade the interest in many body pertur¬ 
bation theory has risen significantly. This is to some 
extent related to the enormous increase in the available 
computer performance, but it is also driven by the real¬ 
ization that many of the presently available density func¬ 
tionals have limited predictive accuracy. Improving den¬ 
sity functionals is a very active field of research in itself. 
In fact, tremendous progress has been made for semicon¬ 
ductor and insulator modelling by the inclusion of exact 
exchangeji^ as well as for weakly bonded systems by in¬ 
cluding either atom centered dispersion corrections^— or 
non-local van der Waals corrections regarding the density 
at two points in space— 2r— However, a unified compre¬ 
hensive, accurate and predictive framework for metallic, 
covalently as well as dispersion driven interactions is hard 
to attain using present density functionals: most avail¬ 
able density functionals require careful evaluation against 
more accurate methods before one can trust them to pre¬ 
dict accurate numbers for a specific material. 

In the quantum chemistry community, such a concise 
hierarchy of methods for evaluating and benchmarking 
more approximate methods, such as density functional 
theory, is well established. The highest rung of this hi¬ 
erarchy is made up by the full configuration interaction 
method, followed by a variety of more approximate meth¬ 
ods. For instance, if the material is dominated by a sin¬ 
gle Slater determinant, the methods of choice are coupled 
cluster methods^ as well as Mpller-Plesset perturbation 
theoryi£ for large band gap systems. Only recently these 


methods have become available for solids —- — 

For solids, calculations using coupled cluster methods 
or full configuration interaction methods are, however, 
exceedingly demanding approaching several 100.000 CPU 
hours for a single material with a few atoms in the unit 
cell. The approach taken in solid state systems is there¬ 
fore often “bottom up”, i. e. starting with an approxi¬ 
mate scheme such as density functional theory and im¬ 
proving upon the description until results compare rea¬ 
sonably well with experiment. The random phase ap¬ 
proximation (RPA) is one promising approach to achieve 
this goal— - In fact, the RPA yields a balanced descrip¬ 
tion of most bonding types, including metallic bonding, 
covalent, ionic, as well as van der Waals (vdW) bonding. 
Initial applications were limited to small molecules* 21 
Although first applications to bulk materials were not 
encouraging,— for many materials results nowadays sur¬ 
pass those from semi-local functionals—The studies 
now span a wide range of applications, including molec¬ 
ular reactions^! rare gas solids^! properties of cova¬ 
lent, metallic and ionic solidsj22r— 2 dispersion forces in 
graphite and between graphene and surfaces layered 
compounds^ adsorption of molecules on surfaces *21 bulk 
ice properties) 33 and many more applications are emerg¬ 
ing. Recent advances include highly efficient implementa¬ 
tions scaling only cubically with system size and linearly 
in the number of k-points^2’ 40 as well as implementations 
scalable to massively parallel computers— 

There is no denying that RPA is not perfect. Among 
all the possible many body diagrams, direct RPA ex¬ 
clusively sums the bubble diagrams. Attempts to in¬ 
clude other kinds of diagrams, for instance higher or- 
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der exchange interactions ; 42 ’ 43 the contribution of single 
excitations^ 4-48 or approximate exchange and correlation 
contributions inspired by density functional theory^ are 
currently vigorously explored research directions. Also 
better starting points than standard density function¬ 
als are explored . 50 ’ 51 Likewise, forces are yet only im¬ 
plemented in two molecular codes ; 52 ’ 53 and they are not 
available in solid state codes. 

The present paper mainly concentrates on the already 
mentioned singles, contributions that arise from determi¬ 
nants where one electron is moved from a state occupied 
in the initial Slater determinant to an originally unoccu¬ 
pied state. If the initial determinant is the Hartree-Fock 
determinant, then the singles are zero in lowest order 
perturbation theory because of Brillouin’s theorem.— If 
one sets out from the Kohn-Sham determinant, the sin¬ 
gles however contribute even in lowest order perturba¬ 
tion theory. In the adiabatic-connection framework this 
contribution was first derived by Gorling and Levy.— In 
the present work, for reasons of consistency we will re¬ 
derive the same term, albeit doing so using a concise 
Green’s function notation, making the derivation alge¬ 
braically simpler (see Sec. Ill Cl) . The main point of our 
work is, however, that we give up the assumption that 
the charge density is kept constant along the coupling 
constant integration (see Sec. im an approximation 
made in most derivations of the adiabatic-connection 
fluctuation-dissipation theorem (AC-FDT) — This ap¬ 
proximation was originally adopted, since density func¬ 
tionals were considered to be very accurate in predicting 
the groundstate density, and because AC-FDT was used 
as a theoretical pathway to derive improved density func¬ 
tionals. However, since AC-FDT is now a computational 
framework, and since density functionals do not always 
yield accurate densities, we feel that this point urgently 
needs to be revised. We are aware of at least two papers 
where this assumption has been given up as well, albeit in 
the first one the point was only made in passing^ and the 
second one only considers small prototypical molecules— 

Our present formal derivation (see Sec. Ill Dl) yields 
results similar to the singles and re-normalized singles 
proposed by Ren et a h 44 ’ 45 ’ 48 Ren and coworkers, how¬ 
ever, formally gave up the adiabatic-connection frame¬ 
work and used standard second order perturbation the¬ 
ory to motivate their singles and re-normalized singles. 
We feel that a concise perturbational framework (here 
the adiabatic-connection framework) is helpful to better 
understand the underlying approximations. We will dis¬ 
cuss that the singles account for the changes of the mean 
field density matrix along the coupling constant integral. 
This is exactly analogous to coupled cluster theory or, 
since coupled cluster theory is just a re-summation of 
certain Goldstone diagrams, standard many body pertur¬ 
bation theory. This insight explains why the inclusion of 
singles increases the bonding between weakly interacting 
fragments, i.e. atoms and molecules. When the singles 
are included the charge density contracts compared to 
the original density functional, and this results in a de¬ 


crease of the Pauli repulsion. We demonstrate this effect 
here for rare gas solids, as well as, the cubic phase of 
ice, the benzene crystal and water adsorption on NaCl. 
For these systems, the predicted cohesive energies and 
lattice constants are, after inclusion of the singles, in 
very good agreement with experiment (see Sec. IIV Al 
IIVBI IIV Cl IIV DD . We also present results for covalently 
bonded as well as metallic systems. Here no improve¬ 
ments are discernible, or rather, the corrections of the 
singles are mostly tiny and do not worsen the already 
excellent agreement with experiment (see Sec. IIV El) . 

II. THEORY 

A. Adiabatic switching with density change 

In this section, we give a brief derivation of the adia¬ 
batic switching, where the correlation energy is obtained 
by switching from the single determinant reference sys¬ 
tem to the many body system of interest. We consider 
the usual adiabatic connection, where the density is kept 
constant along the pathway, as well as adiabatic switch¬ 
ing allowing for a change of the density. 

In the adiabatic-connection framework, it is common 
practice to switch from the purely local, in real space 
multiplicative potential V\ to the exact many electron 
Hamiltonian H\— 

H\ =^/r(r i ) + ^^u(r l ,r') + ^V A (r i ). (1) 

* &3 * 

Here v is a two-electron operator, typically the Coulomb 
kernel l/|r — r'|, and h is an arbitrary one electron op¬ 
erator, typically the kinetic energy operator A and the 
potential of the ionic cores or some other external one- 
electron potential V ext 

k r) = -iA r + V ext (r), 

where we have assumed atomic (Hartree) units. The cor¬ 
relation energy £’ lot is defined as the difference between a 
single Slater determinant To evaluated using the orbitals 
at zero coupling (A = 0) and the exact many-electron 
wave function Ti evaluated at full coupling (A = 1): 

E%* = (* 1 \H 1 \* 1 )-(*o\Hi\*o)- (2) 

At zero coupling, the Hamiltonian is chosen to be the 
usual Kohn-Sham Hamiltonian, 

(“ A + V ext + V 0 KS ) |0„) = e n \<t> n ), (3) 

where V 0 KS i s the standard Kohn-Sham potential (we 
keep the subscript 0 to indicate zero coupling). Note 
that the Kohn-Sham potential includes exchange, corre¬ 
lation, and the Hartree contributions. In the following, 
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we will also use the shorthand | n) = \(f > n ), and the indices 
always refer to occupied one-electron orbitals, 
whereas the indices a, b, c,... refer to virtual (unoccupied) 
one-electron orbitals of the Kohn-Sham Hamiltonian. 

Subtracting and adding K 0 KS to the second term on the 
r.h.s of 0 one can rewrite the correlation energy as 


Linear switching is also exact (as long as no phase tran¬ 
sition is encountered) and yields results identical to stan¬ 
dard Rayleigh-Schrodinger perturbation theory where 
the perturbation is also switched on linearly. 

Linear switching yields for the second line in Equ. © 
a simple correction to the correlation energy 


El ot = ('F 1 |iL 1 |'F 1 ) - (^o|iL 0 |^o) 

- i(tf 0 |fi|*o) + <^o|V" 0 KS |^o> 

E™ = x\^\^x)dX W 

- J(4-o|^|^o) + («'o|V r 0 KS |«'o)- 

In going from the first to the second equation, the 
“Hellman-Feynman” theorem has been used, i.e. it is 
assumed that is the groundstate wave function of 
the Hamiltonian H\, and terms involving the derivative 
of the wave function d'&x/dX thus vanish. The first 
term in the second line is exactly the Hartree and ex¬ 
change energy evaluated for Kohn-Sham orbitals. Insert¬ 
ing the derivative of the Hamiltonian dH\/d A immedi¬ 
ately yields (compare Equ. dTJ ) 

K ot = l /Va|«|*a )d\ - J<*oM*o> + 

,. (5) 

J ^A|^|vhA)dA + (T 0 |K 0 KS |^o). 

We note that a similar term is given in the appendix of 
Ref. [56j without a clear derivation. The most common 
way to realize the coupling integral is to chose the poten¬ 
tial V\ such that the density remains constant along the 
coupling constant integral. 20 Then the expectation value 
of the density operator yields the (constant) groundstate 
density (’F A |h(r)|'I'A) = n(r), and the first term in the 
second line can be integrated yielding 

= <*o|Vi - K KS |T 0 >, 

canceling the last term in the second line of Equ. ©• 
Since the additive potential V\ must be zero at full cou¬ 
pling V\ = 0, one obtains the standard equation for the 
AC-FDT: 

Ec = lJ a - ^<*o|«|* 0 >, (6) 

which is only valid if the density is constant along the 
integration pathway , whereas the full version is obviously 
given by Equ. ©• 

To obtain a simplified version of the full equation, one 
can switch off the Kohn-Sham potential linearly i.e. 

V x = (1 - A)K 0 ks . (7) 

Linear switching was first considered by Harris and 
Jones^ and subsequently discussed in Refs. |H and H3- 


A E c = 



d\ J d 3 rn A (r)b 0 KS (r) + J d 3 r n 0 (r)K 0 KS (r) 
dX J d 3 r (n A (r) - n 0 (r)) K 0 RS (r). 

( 8 ) 


Here nx is the charge density at coupling A. The total 
correlation energy is then given by the sum of Equ. © 
and Equ. ©. 


B. Fluctuation-Dissipation theorem 


It is common to rephrase the correlation energy in 
Equ. © using the fluctuation-dissipation theorem. The 
derivation is sketched in appendix A, however, similar 
derivations can be found elsewherei 20 ’ 56 ’ 57 ’ 59 The final 
result becomes: 


E c = 


([n\vnx] - [riowio]) dX 


If 1 

2 J Ti '[(xa(0 ) + n A -xo(0 )-n 0 )v]dX. 


(9) 


Here we have introduced short-hands 


[nvn] = J n(r)v(r, r')n(r')d 3 rd 3 r' 

T.-[(x + n)v] = J M r.r’) + »(r)J(r - r^.r^ri 


( 10 ) 


In the equations above, xx is th e reducible polarizability 
(or density fluctuation response function) of the many 
electron system at coupling A: 

XA(T,r,r') = -(ifx\T[6n(T,r)6h(0,r')]\ipx)- (H) 


T is the time ordering operator, and 8h the density fluc¬ 
tuation operator Sri = h — (h)x- The polarizabilities 
are evaluated at a small negative infinitesimal imaginary 
time r — } 0~. Note that throughout the paper, r is 
the imaginary part of an imaginary time t = — ir. The 
polarizability in imaginary time can be related to the po¬ 
larizability in imaginary frequency ico by a Fourier trans¬ 
formation: 


Xa(t) 

Xa(*w) 


h 


— oo 
poo 


du ex.p(—iuJT)x\(iaj) 
dr exp(iwr)x A (r). 


( 12 ) 
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At zero coupling, Xo becomes the independent particle 
polarizability Pq of the Kohn-Sham system, and can be 
written in terms of the one particle Green’s functio n 56 ’ 60 
We define the independent particle polarizability at cou¬ 
pling A generally as 

P\(t, r,r') = G\(t, r, r')G\(—r, r / , r), (13) 


where G\ is the one-particle Green’s function of the A 
interacting system H\. At zero coupling, i.e. for the 
Kohn-Sham case, the Green’s function in imaginary time 
Go are given by (/i is the chemical potential of the elec¬ 
trons, i are occupied and a are unoccupied one electron 
orbitals): 


Go(r) 


£il*X*l exp(-(ei - h)t) r<0 , , 

X) 0 |a)(a| ex p(-(ea-M)'r) ^>0 


Before continuing, we note that at small A and neglect¬ 
ing fluctuations, the term in the integral on the second 
line in Equ. © can be also written as an integral of the 
exchange energy 

l 

7 *(r, r') 7 A (r , ,r){)(r , ,r)fiA, (15) 

where 7 ^(r,r') is the one-particle density matrix at cou¬ 
pling A. This follows from expanding the S function in 
one-electron orbitals <5(r,r / ) = J]mgaii( r l TO )( TO l r, )i an d 
inserting the independent particle approximation for the 
polarizability (compare Eaus.fTdlandfTdll. Hence, in Equ. 
© the first line represents the change of the Hartree en¬ 
ergy along the coupling constant integration, whereas the 
second line accounts for the change of the exchange en¬ 
ergy (for uncorrelated wave functions). 

Returning now to Equ. m, we note that replacing 
X\ by P\ in Equ. © neglects important many body 
effects related to changes in the Hartree or exchange po¬ 
tential. In many body perturbation theory, these can be 
included exactly by solving the Bethe-Salpetlier equa¬ 
tion for the polarization propagator—f~ 64 The simplest 
approximation to the BSE equation is the common RPA 
approximation 



Xa(w) = PxM + P\(u)\vx\{u), (16) 


which only includes the Hartree-related ring or “bubble” 
diagrams. In time dependent DFT, the related equation 
is 


Xa(w) = PoM + Po(w)( Xv + / xc ( A))xa(w), (17) 

where / xc (A) accounts for all correlation effects, includ¬ 
ing the change of the independent particle polarizability 
along the coupling path. 

If the density is kept constant (n\ = no) and one ap¬ 
proximates P\(ui) = P 0 (u;) in Equ. (fThl) . one obtains for 
Equ. © the direct RPA correlation energy^ 

1 r°° 

E c PA= 2 Wo dwTr Hl--PoM^)+PoM^]. (18) 


Here and in the following, we often drop position indices 
and implicitly assume integration over the missing spatial 
coordinates: 

(AB)( r,r") = J A(r,r')B(r',r")d 3 r'. 


C. Singles contribution in Gorling-Levy 
perturbation theory 


Before considering density changes along the coupling 
constant integral, we will briefly derive the singles expres¬ 
sion in the AC-FDT for a fixed charge density in order 
to compare with this equation later. 

In standard RPA, one neglects that the independent 
particle polarizability changes along the coupling con¬ 
stant integration, as the one-particle Green’s function 
changes as A changes (see Equ. (fl9l) below). Oddly this 
term is not often considered in the AC-FDT framework, 
although Gorling and Levy have already highlighted its 
relevance (albeit not in a Green’s function formalism). A 
partially analogous discussion using the Green’s function 
formalism can be found in e.g. Ref. |H andl66l. 

At small coupling A, the change of the one-particle 
Green’s function is described by the first order term in 
the Dyson equation (see e.g. Ref. [67l ): 

/ OO 

dr Go(t)(AV ks + V x )G\{—r) 

-oo 

/ OO 

dr G 0 (t)(AE ks + V x )Gq(—t), 

-oo 

(19) 

where Ab KS (r)i(r,r') is the change of the local Kohn- 
Sham potential, and V x (r,r') is the exact non-local ex¬ 
change potential evaluated using the DFT orbitals at cou¬ 
pling constant zero: this term originates from switching 
on the exact many body potential, which in lowest order 
is equivalent to switching on the non-local exchange. If 
one performs the coupling constant integral keeping the 
density fixed, then the change of the Kohn-Sham poten¬ 
tial must be chosen such that the density remains exactly 
constant to the original density at A = 0. This require¬ 
ment can be written as: 


n\(r) — n 0 (r) = lim (G a (t, r, r) - G 0 (r, r, r)) = 0, 

T—>0 

( 20 ) 

as the diagonal of the Green’s function at r —»• 0 
is just the charge density. From Eq. C© one obtains 
the linearized Sham-Schliiter equation for the potential 

AI/KS^ 


dr G 0 (r){AV Kb + V x )G 0 {-r) = 0. 


( 21 ) 


Since the Kohn-Sham potential is local, one can factor 
out Pq(uj = 0) = / g?tGo(t)Go(—t) (see Equ. d©) and 


5 




^7 


FIG. 1. Goldstone diagrams corresponding to the singles. 
The sign is given by the number of closed Fermionic loops 
(negative for first and last diagram, positive for the two in the 
middle). The exact exchange potential is density dependent 
(bubble) and indicated by a broken line. 


obtain: 

/ OO 

cItGo(t)V x Gq(—t) . 

-oo 

( 22 ) 

This is just the local exact-exchange optimized-effective 
potential (EXX-OEP) P EXX . In summary, if the den¬ 
sity is supposed to remain constant along the coupling 
constant, the first order change of the local Kohn-Sham 
potential is exactly given by the exact-exchange OEP po¬ 
tential. This is not surprising, since the defining property 
of the Sham-Schliiter equation is that the density from 
the non-local exchange potential and the effective local 
potential must equal each other. 

One can now calculate the change of the correlation 
energy by inserting the first order change of the Green’s 
function into the expression for the correlation energy 
in the AC-FDT. Approximating x\ ~* P\ = G\G\, in¬ 
serting m and £0 into the second line of Equ. m, 
and rewriting the S function as a sum over states, or 
alternatively, starting from Equ. m and identifying 
7a = Ga(0 - ) yields a change of the correlation energy of 
(yet omitting the integration over A): 

~ l°° (tt[G 0 (t)(V x - P exx )G o (-t)Go(0-)D 

+ Tr[G 0 (0-)Go(T)(l7 - P exx )G o (-t)0]) dr. 

(23) 


Because of the trace, the second term yields just the 
complex conjugate of the first term. Furthermore, the 
integral over the positive half-plane of r gives the same 
value as that over the negative half-plane (r < 0). Fi¬ 
nally, the term — Go(0~)0 corresponds to the non-local 
exchange potential V x . Inserting the defining equation 
for the one-particle Green’s function (1141) and perform¬ 
ing the integration over r and A yields what is commonly 
referred to as the singles contribution: 


E? = - X 

aGvirt, 


l(i|t> a: -t> EXX |a)| 2 
e a ^ G 


(24) 


Note that one needs to add Equ. (l^Tl) to Eq. (E51) mul¬ 
tiplied by P EXX (r) and integrated over r to derive this 


convenient equation. The corresponding (time-ordered) 
Goldstone diagrams are also shown schematically in Fig. 
ID A few comments are in place here. The term has been 
first derived by Gorling and Levy^ In their derivation 
it is, however, not obvious that this term describes the 
changes of the correlation energy from the one-particle 
density matrix (which equals the one particle Green’s 
function at r = 0“) along the coupling constant integral. 
In fact, most standard AC-FDT calculations neglect the 
term. Only the RPAX (RPA including eXchange) follow¬ 
ing Gorling and coworkers accounts for this term (often 
referred to as EXX-RPA)^ In these formulations, how¬ 
ever, the change of the Green’s function is accounted for 
by recasting it (as well as the particle-hole ladder dia¬ 
grams) into an effective exchange kernel f xc ( A) ss A f x 
for the polarizability (compare Equ. [T71) j 66 i 69 

We conclude, the singles account for the change of the 
density matrix along the coupling constant integral. For a 
constant density, obviously only changes of the exchange 
and kinetic energy can be included. In the next step, we 
will also allow for changes of the charge density along the 
coupling constant integral. 


D. Singles contribution with density changes 

We now derive the singles contribution to the corre¬ 
lation energy for the case when the density does not 
stay constant during the coupling constant integration. 
Changes of the charge density are most easily accounted 
for by linearly switching off the Kohn-Sham potential (see 
Equ. 0) and linearly switching on the exact many body 
interaction. In principle, this makes the derivation even 
less involved, since the determination of the local exact 
exchange potential is no longer required. We first again 
derive the expression in lowest order, where the Green’s 
function is now given by 

/ OO 

dr'G 0 (tO(P HF - P 0 KS )Go(-t'). 

-oo 

(25) 

In the lowest order, the change of the potential is now 
the difference between the Hartree-Fock potential P HF = 
V H + V X , the sum of the Hartree and exchange potential, 
and the original Kohn-Sham potential, which is adiabat- 
ically switched off. As before, the change of the density 
matrix is given by the change of the Green’s function 
G(t — > 0“) (compare Equ. (f20|) 4 . In this case, the first 
and second line of Equ. 0 can be combined to yield: 

A r°° 

- / drTr[G 0 (r)(I> HF -I/ 0 KS )Go(-r)K HF ]+c.c. (26) 

The first line of Equ. 0 yields the Hartree-potential 
times the change of the density, whereas the second line 
yields the exchange potential times the change in the 
density matrix, in sum the change of the density matrix 
times P HF . 
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After performing the coupling constant integral this 
yields 


E c 


y, (z|K HF -K 0 KS |a)(q|F HF |i) 
e a - e» 

iGvirt, 
i(z occ 


(27) 


This term needs to be combined with the term A E c 
given in Equ. m, which can be calculated by inserting 
Equ. (fTfll) and performing the r and A integration. Both 
contributions combined yield a simple term 


term Go(0 ) is integrated over so that we obtain 

r i 


Tr[G A (0-)(/i HF -h KS )\d\ 
-Tr[G 0 (0-)(h HF -h KS )]. 


(31) 


As the one-electron Green’s function G\ (0 ) is the exact 
Green’s function of the one-electron Hamilton operator 
hx = h KS + A(/jHF _ h KS ), one can rewrite the first line 
as (Hellman-Feynman theorem): 


E c se + A E c 


y |(z|E HF -E 0 KS |a )| 2 

iGvirt, f ( ' 

iG occ 


(28) 


which exactly corresponds to the singles suggested by 
Ren et ah 44 ’ 45 Here, we have performed the deriva¬ 
tion concisely within the AC-FDT framework instead 
of Rayleigh-Sclrrodinger perturbation theory, and, as it 
must be, both are exactly equivalent. As in the previ¬ 
ous paragraph, the singles account for the change of the 
mean field exchange energy. However, now they also in¬ 
clude the change of the mean field Hartree energy along 
the coupling constant integration. Here and in the fol¬ 
lowing, we define the mean field as the contributions that 
stem from the one-particle Green’s function and the re¬ 
lated density matrix. 

In the renormalized singles of Ren and coworkers^ also 
higher order contributions are accounted for. However, 
it is not entirely straightforward to generalize our results 
to include higher order terms in A, and to continue, we 
make one crucial approximation. Let us introduce this 
approximation for the density term, which can be written 
as (compare Equ. ©) 

r {[n x vn\] - [n Q vno]) = 

\ (29) 

~[(n\ - n 0 )v(n\ +n 0 )] « [(n A - n 0 )vn 0 }. 

In the second line, we assume that the differences be¬ 
tween n A and no are small so that we can approximate 
the sum by 2no- Analogous manipulation is possible 
for the term involving the polarizability, if the full po¬ 
larizability is approximated by the independent particle 
approximation Equ. 1H3I) (x A —t P\ — G\G\). After 
collecting all terms, adding Equ. <[8|), and noting that 
n(r) = G(0 _ ,r, r), one obtains the approximate renor¬ 
malized singles correlation energy 

~ [ Tr[(G A (0-)-G 0 (0-))(t> HF -H o KS )]dA. (30) 

Jo 

It is fairly straightforward to backtrack this into a sim¬ 
ple total energy relation (essentially reversing the steps 
introduced in Sec. Ill Al) . We first note that the difference 
between the Hartree-Fock potential and the Kohn-Sham 
potential equals the difference in the corresponding one- 
electron Hamiltonians h HF — h KS . Next, the constant 





d\ = Tr 


Gi(0 _ )h 1 


G o (0 )ho 


Combining this with the second line in Equ. (EHl) {ho = 
h KS , hx = h HF ) yields 

E r c SE = Tt[7hf^hf — TdftAhf]- (32) 

Here 7 hf = Ghf(0 _ ) is the Hartree-Fock density matrix, 
determined for the Hartree-Fock Hamiltonian /ihf , where 
the Hartree-Fock potential is set up with DFT-orbitals. 
This is exactly the “single-shot” Hartree-Fock energy: it 
can be calculated by diagonalizing the HF Hamiltonian 
(set up with DFT-orbitals), summing the eigenvalues of 
the occupied states and subtracting the original diagonal 
part of the HF Hamiltonian for the occupied manifold 
evaluated using the original DFT orbitals (i DFT ): 


Ef E = £ e? F - (i DFT |h HF |* DFT ). (33) 

i 

This prescription is one central result of the present pa¬ 
per. 

In all practical calculations, we found the value of the 
single shot HF energy to be exceedingly close to the 
renormalized singles introduced by Ren and coworkers4^ 
For diamond the difference is on the order of 0.3 meV 
(with the singles being of the order of 0.3 eV). Even for 
small band gap systems, such as metallic Pd or Al, differ¬ 
ences hardly ever exceed 1 % and are entirely irrelevant 
when evaluating relative energies or lattice constants. 
The relation of our equation to the singles of Ren et al. is 
fairly straightforward to see. Ren et al. essentially renor¬ 
malizes the propagator G in the occupied-occupied block 
as well as the unoccupied-unoccupied block, by diago- 
nalization of these sub-blocks using the HF Hamiltonian 
/ihf- This step is crucial, since the one-electron eigen¬ 
values are renormalized to the HF-eigenvalues; if LDA 
eigenvalues were used in the evaluation of the singles in 
Equ. (1281) , the response of the system to the change of the 
potential from Kohn-Sham to HF would be strongly over¬ 
estimated. Ren then calculates the change of the mean 
field energy in second order with the DFT eigenvalues in 
Equ. (l28l) replaced by the renormalized HF eigenvalues. 
As opposed to this prescription, we also “renormalize” 
the propagator in the occupied-unoccupied block. With 
the present derivation at hand, there is no obvious reason 
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why not to chose the simpler prescription of the present 
work. Since results obtained using the renormalized sin¬ 
gles (rSE) are in practice indistinguishable from results 
obtained using the single-shot HF energy change, all cal¬ 
culations here use the single-shot HF energy change, but 
are nevertheless labeled as “rSE”. 


E. Singles contribution beyond the Hartree-Fock 
description 

An obvious extension to the renormalized singles ap¬ 
proach is to use the full RPA density matrix instead of 
the HF density matrix to estimate the change of the mean 
field energy: 

-EGWse = Ti'[7rpa/ihf — 7dft^hf], (34) 

where the 7 rpa is the RPA density matrix. In the singles 
derived in the previous section, it was assumed that the 
density matrix of the interacting system is well approx¬ 
imated by the Hartree-Fock case, which seems a crude 
approximation. Given the excellent performance of the 
RPA for total energies and band gaps, the RPA density 
matrix, however, should approximate the density matrix 
of the real interacting system very accurately. At this 
point, we have, however, no entirely concise derivation 
for the term, although the physical interpretation is clear: 
it should account approximately for the change of the 
mean-field kinetic, Hartree and exchange energy along 
the coupling constant integral. 

To evaluate the RPA density matrix, we calculate the 
RPA Green’s function 

Grpa(*w) = G 0 {iu) + G R p A (*w)(E(iw) - V* s )G 0 {iu;). 

(35) 

Here E is the self-energy in the GW or random phase 
approximation (the two approximations are synonymous) 

S(t, r, r') = —Go(t, r, r')W (r, r', r), (36) 

and Go and W are the Kohn-Slram Green’s functions 
(fUf!) and the screened potential evaluated using Kohn- 
Sham polarizabilities^ V^ s j s ^he exchange-correlation 
contribution to the Kohn-Slram potential. 

To evaluate the RPA density matrix numerically, we 
determine the interacting Green’s function at full cou¬ 
pling Grpa(*w) using Equ. (l35l) . transform it to the imag¬ 
inary time at a small negative infinitesimal r —> 0~ 
to obtain the one-particle density matrix for the RPA, 
7 RPA(r,r / ) = Grpa( 0~, r, r'), and finally diagonalize the 
density matrix to obtain the natural orbitals 

7RPA = ^2 \ m ) fm { m \. 
m 

In the present work, we use one more crucial approxi¬ 
mation: instead of using the occupancies of the actual 
density matrix, we keep the original occupancies evalu¬ 
ated on the level of density functional theory (1 and 0 


for insulators). This has two reasons: first the density 
matrix evaluated from the Green’s function (1351) is not 
particle conserving, be. the number of electrons deviates 
from the original number.® 7 Only if the Green’s function 
in (1351) and (1361) is iterated to self-consistency, the par¬ 
ticle number is conserved. To test our present code, we 
have also iterated the Green’s function in both equations 
to selfconsistency, and in that case, the electron number 
is indeed exactly conserved. However, such calculations 
are fairly expensive and difficult to apply routinely. They 
would also most likely require us to combine it with a 
different treatment of the fluctuation terms beyond the 
standard RPA treatment as used here* 67 

The second reason is based on the definition of sin¬ 
gles in quantum chemistry. A density matrix with oc¬ 
cupancies 0 and 1 corresponds to a single Slater deter¬ 
minant. We, therefore, approximate the density matrix 
by the “best” single Slater determinant approximating 
the correlated RPA density matrix. This is consistent 
with Brueckner coupled cluster orbitals, 70 which are ob¬ 
tained by performing a rotation between the occupied 
and unoccupied manifold to determine an optimal refer¬ 
ence single Slater determinant. The rotation is chosen 
to remove all “singles” contributions from the correla¬ 
tion energy. In quantum chemistry, the fluctuations that 
cause the fractional occupancies in the density matrix 
are, in fact, not included in the singles. For instance, 
in the coupled cluster theory, fluctuations are accounted 
for by double excitation operators (doubles). Per defini¬ 
tion, singles are diagrams ending in a single excited de¬ 
terminant with one hole in an orbital i and an additional 
electron in an orbital a (compare Fig. [[]). We essentially 
follow the quantum chemistry convention in partitioning 
the correlation energy into a “singles” term and fluctu¬ 
ation terms described entirely by Equ. (fl8l) . We finally 
note that the RPA orbitals constructed in this manner 
are similar but not identical with the Brueckner RPA di¬ 
rect ring coupled-cluster orbitals suggested by Moussa^ 


III. COMPUTATIONAL SETUP 

The present calculations were performed using the 
VASP code. 7 - We used the new implementation of the 
RPA routines as documented in Refs. IH and The 
code has been extended to allow for (self-consistent) 
calculations of the one particle Green’s function in the 
random phase approximation i.e. solving Equ. m 
and m- Here only single shot calculations are per¬ 
formed, by calculating the self-energy once at full cou¬ 
pling and determining the RPA Green’s function and 
RPA density matrix (Equ. (155l) l. The singles contribu¬ 
tions are evaluated either using the single shot HF den¬ 
sity matrix (Equ. (13^1 ) ) or the single shot RPA density 
matrix (Equ. (l34l) f. We term the two results either rSE 
(renormalized singles) or GWSE, respectively. For the 
fluctuation term we use the standard random phase ap¬ 
proximation as defined in Equ. USD- Furthermore, we 


TABLE I. PAW Potentials and plane wave cutoffs for the or¬ 
bital basis used for the individual calculations (as distributed 
with vasp.5.4.1). The shorthand, _sv indicates that the up¬ 
per most core shell was included as “valence” orbitals (see 
explanation in Sec. II V El) . For the norm-conserving case, 
the enumerated potentials were replaced by the same poten¬ 
tial with an appendix _nc, and for N, O, and F the small 
core potentials N_h_GW, 0_h_GW, and F_h_GW were used. 
These potentials are almost norm-conserving. For the norm- 
conserving potentials always the default cutoff was used. 



applied potentials 

cutoff (eV) 

Ne 

Ne_GW 

1000 

Ar 

Ar_GW 

600 

Kr 

Kr_GW 

500 

ice 

O.GW H.GW 

600 

benzene 

C.GW H.GW 

600 

water on NaCl Najsv.GW Cl.GW O.GW H.GW 

400 

Na 

Na_sv_GW 

350 

A1 

Aljsv.GW 

534 

Rh 

Rh_sv_GW 

351 

Pd 

Pdjsv.GW 

356 

Cu 

Cu_sv_GW 

509 

Ag 

Ag_sv_GW 

460 

C 

C_GW_new 

414 

Si 

SLGW 

245 

Ge 

Ge_sv_GW 

533 

LiF 

Li_sv_GW F.GW 

487 

LiCl 

Li_sv_GW C1_GW 

433 

NaF 

Na_sv_GW F.GW 

487 

NaCl 

Na_sv_GW CLGW 

341 

MgO 

Mg_sv_GW O.GW 

434 

SiC 

SLGW C_GW_new 

414 

AIN 

Aljsv.GW N_GW_new 

547 

A1P 

Al_sv_GW P.GW 

534 

AlAs 

Al_sv_GW As_sv_GW 

539 

GaN 

Ga_sv_GW N.GWjiew 

420 

GaP 

Ga_sv_GW P.GW 

404 

GaAs 

Ga_sv_GW As_sv_GW 

539 

InP 

Injsv.GW P.GW 

476 

Tn As 

In_sv_GW Asjsv.GW 

539 

InSb 

Injsv.GW Sbjsv.GW 

484 


often show the results of the random phase approxima¬ 
tion (evaluated using DFT orbitals) combined with HF 
exchange energies. In this case, the exact exchange en¬ 
ergy evaluated using DFT orbitals is replaced by the self- 
consistent HF energy. The DFT orbitals are always de¬ 
termined using the PBE functional, where PBE stands 
for the usual Perdew-Burke-Ernzerhof functional. 7 - The 
potentials are generally equivalent to the GW potentials 
distributed with vasp .5.4.1. If not otherwise stated the 
default plane wave cutoffs were used for the calculations. 
For the response function, the cutoff was set to 2/3 of 


TABLE II. Lattice constants (in A) as well as cohesive ener¬ 
gies (in meV) of rare gas solids calculated using the RPA and 
different variants for the kinetic, exchange and Hartree energy 
contributions. The experimental data have been corrected for 
zero point vibration energies as determined by accurate quan¬ 
tum chemical methods.— 


volumes 

EXX+RPA 

+rSE 

+GWSE 

HF+RPA 

EXP 

Ne 

4.35 

4.32 

4.38 

4.28 

4.30 

Ar 

5.33 

5.24 

5.25 

5.20 

5.25 

Kr 

5.68 

5.61 

5.61 

5.56 

5.63 

energies 

EXX+RPA 

+rSE 

+GWSE 

HF+RPA 

EXP 

Ne 

19 

35 

30 

36 

26.2 

Ar 

69 

88 

87 

93 

87.9 

Kr 

97 

119 

119 

126 

121.8 


the cutoff for the orbitals. If the energy-volume curves 
were non smooth using 4x4x4 k-points, the plane wave 
energy cutoffs were increased by 30 %, roughly doubling 
the number of included plane waves. The applied compu¬ 
tational setup and potentials are summarized in Tab. [T] 
To determine the equilibrium volumes, the volumes 
were typically varied in steps of 5% by ±15% around the 
experimental equilibrium volume, and the Murnaghan 
equation of state was fitted to the data points— 


IV. RESULTS 
A. Rare gas solids 

Rare gas solids constitute a prototypical test case for 
van der Waals bonded solid state systems. Although, 
the standard RPA performs reasonably well for rare gas 
solids,— one observes that the binding energy is quite sig¬ 
nificantly underestimated, in particular for He. Inclusion 
of singles has remedied this issue for the rare gas dimers, 
and one would expect that this also applies to solids—^ 

We first note that quantum chemical coupled cluster 
calculations at the level of CCSD(T) (coupled cluster 
with singles and doubles and perturbational triples) us¬ 
ing the method of increments and including up to four 
body interactions yield essentially exact results within a 
few /rHartree.' 75 ' 76 To compare with experiment we have 
corrected the experimental data for zero point vibration 
energies in both the cohesive energy, as well as in the 
lattice constants— 

The present calculations have been performed using 
relatively high plane wave energy cutoffs of 1000 eV, 
600 eV and 500 eV for Ne, Ar and Kr in order to ob¬ 
tain smooth energy-volume curves. 4x4x4 k-points were 
used, although already 3x3x3 k-points yield practically 
identical results. The PAW potentials are approximately 
norm-conserving to avoid errors in the vdW contributions 
from excitations into high lying unoccupied orbitals. 77 As 
already mentioned, standard RPA combined with exact 
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exchange yields at best modest agreement with experi¬ 
ment (see Tab. ED- Specifically, the equilibrium volumes 
are overestimated and the binding energies are, as al¬ 
ready noted above, about 20 % too small, with the errors 
being particularly large for Ne. Inclusion of the singles 
contributions from HF (rSE) and using the RPA density 
matrix (GWSE) yields a clearly improved agreement with 
experiment in particular for Ar and Kr. Ne is less sat¬ 
isfactory. For Ne, the binding energy significantly over¬ 
shoots if the singles are evaluated on the Hartree-Fock 
level (rSE). This improves, if the RPA density matrix is 
used, however, including GW-singles worsens the volume 
even compared to RPA. We believe that this is a result of 
the PBE approximation being particularly inaccurate for 
large band gap system such as Ne. This is for instance 
exemplified by the fact that the binding energy almost 
doubles, when the exact exchange evaluated using DFT 
orbitals is replaced by the exact exchange evaluated us¬ 
ing Hartree-Fock orbitals (HF+RPA). For Ar and Kr the 
changes are typically only 25 %. Also other metrics indi¬ 
cate that the error in the orbitals is significant using the 
PBE functional for this case. Hence, single shot correc¬ 
tions, be it RPA+rSE or RPA+GWSE, yield less reliable 
results then for Ar and Kr. For Ar and Kr, the per¬ 
formance of RPA+rSE and RPA+GWSE is remarkably 
good, approaching that of state of the art quantum chem¬ 
istry methods. Furthermore, the differences between rSE 
and GWSE are small (except for Ne), which is related to 
the fact that all these systems screen relatively weakly. 
Therefore, the RPA density matrix is very close to the 
HF density matrix (see below). 


B. Ice 

Ice is another system, where singles are expected to 
have a significant impact. As for rare gas solids, the 
PBE charge density is too spread out and replacing the 
exact exchange evaluated using PBE orbitals by the ex¬ 
act Hartree-Fock exchange increases the binding energy 
by 100 meVi^ The predicted volumes using EXX+RPA 
and HF+RPA almost exactly bracket the experimental 
values. 

Here we only concentrate on one ice phase, the low¬ 
est energy cubic phase of ice, I c (a), with a ferroelectric 
order. In our previous study, we found this phase to be 
practically iso-energetic with the ferroelectrically ordered 
ice XI in the Cmc2i space group. Ice XI is a proton or¬ 
dered variant of the common proton disordered phase of 
hexagonal ice. In the present study, we used PBE re¬ 
laxed structures and identical potentials as in our previ¬ 
ous study^ However, the cutoff was increased from 450 
eV to 600 eV. This generally results in smoother energy- 
volume curves and changes the calculated volumes by 
about 0.5 %: because of the increased cutoff and reduced 
noise in the calculated data, the present calculations are 
slightly more accurate. 

As already observed in our recent work, the 


TABLE III. Equilibrium volume per molecule (in A 3 ) as well 
as cohesive energy (in meV) with respect to an isolated H 2 O 
molecule for cubic ferroelectrically ordered ice I c (a). 



EXX+RPA 

+rSE 

+GWSE 

HF+RPA 

EXP 

volume 

32.94 

31.76 

32.03 

31.38 

32.105“ 

energy 

536 

630 

620 

670 

610 


“ Ref. Zl 


EXX+RPA underbinds, whereas combining the RPA 
with Hartree-Fock energies overestimates the binding en¬ 
ergies (see Tab. IIIII) . Including the singles in the HF 
approximation improves the description already signif¬ 
icantly, although the results are still too close to the 
HF case. In this case, the singles evaluated using the 
RPA-density matrix (GWSE) yield results very close 
to experiment and practically identical to the diffusion 
Monte Carlo data, which gives an equilibrium volume per 
molecule of 31.7 A 3 and a binding energy of 605 meV for 
hexagonal ice^ Since hexagonal ice has a 0.2 A 3 smaller 
volume and a 5 meV reduced binding energy due to dis¬ 
order (compare Table III in Ref. [38D. the results are vir¬ 
tually on top of the DMC data at a tiny fraction of the 
compute cost. 


C. Benzene 

Benzene is the simplest aromatic molecule and there¬ 
fore it represents an interesting reference point for the 
study of more complex molecular solids. The delocalized 
character of 7r bonds in benzene makes an accurate calcu¬ 
lation of interactions or cohesive energies difficult, both, 
in the gas phase, as well as in the condensed phases. Sim¬ 
ple methods, such as MP2, lead to overestimated binding 
energies and more involved schemes need to be used to 
obtain accurate results. For example, for crystalline ben¬ 
zene Wen and Beran computed a post-MP2 correction of 
10.4 kJ/mol (108 meV)^ reducing the MP2 cohesive en¬ 
ergy considerably. Li et al. estimated the RPA cohesive 
energy to be 47 kJ/mol or 487 meVsignificantly below 
the recently revised experimental estimate of —55.3 ± 2.2 
kJ/mol (about 573 meV)42 

We have used the geometry of the crystal and monomer 
optimised with the optB88-vdW functional and per¬ 
formed both the extrapolation to infinite k-point mesh 
and infinite basis sets. To obtain the reference energy of 
the monomer, we also performed calculations at differ¬ 
ent volumes and extrapolated to the infinite cell volume. 
The RPA and EXX calculations used up-to 3x3x3 and 
6x4x6 k-point meshes, respectively. 

The data are summarised in Table HV1 Our EXX+RPA 
value is similar to the previous calculation of Li et al. and 
underestimates the reference by about 110 meV. The ef¬ 
fect of the single excitations is rather small, the bind- 
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TABLE IV. Cohesive energies of benzene crystal obtained 
with different methods and compared to the reference value 
derived from experimental data. 


Method 

E co h (meV) 

EXX+RPA 

-463 

+rSE 

-503 

+GWSE 

-518 

HF+RPA 

-520 

Reference 

-573 ± 23 


ing energy increases by 40 meV for the HF based sin¬ 
gles and by about 55 meV when we use the GW singles. 
This means that the final error is halved compared to 
the original RPA calculation. The hybrid scheme, where 
self-consistent Hartree-Fock is used for the mean field 
part, gives results in a better agreement with the refer¬ 
ence data, underestimating it by about 43 meV. We note 
that such a performance is in agreement with the results 
of Ren et al. 4^ who found that RPA underestimates in¬ 
teraction energies for benzene dimers both in the stacked 
and in the T-shaped geometries. Moreover, including 
single excitations did not improve the results substan¬ 
tially. This clearly shows that singles have limitations. 
We believe that the main error in this case stems from 
an inaccurate description of the delocalized n electrons 
in PBE. 


D. Water on NaCl 

To explore the accuracy of the RPA and the singles 
corrections for adsorption on surfaces, we have calculated 
the adsorption energy of a water molecule on the NaCl 
surface. This is a prototypical system and one for which 
an estimate of the adsorption energy of water has been 
obtained— However, this reference data was calculated 
using an embedded cluster approach to obtain the corre¬ 
lation energy. This makes a direct comparison to our cal¬ 
culations, that are necessarily a finite coverage, difficult. 
Therefore, we opted not to use this reference data and 
instead obtained an estimate from MP2 calculations for a 
small system. To reduce the computational cost of MP2, 
we, furthermore, restrict the study to a small p(lxl) su¬ 
percell with a single water molecule and two layers of 
NaCl and use that consistently for both MP2 and RPA. 
This corresponds to a high coverage, with one molecule 
per one surface sodium atom. Moreover, the reference 
surface and isolated molecule have the same geometry as 
in the adsorbed structure, and the same simulation cell is 
used for all cases to allow for efficient error compensation. 

The interaction energy depends only weakly on the cut¬ 
offs chosen for orbital and auxiliary plane-wave basis sets. 
For MP2 we set the cutoff for the orbital basis and aux¬ 
iliary polarizability basis to 350 eV and 450 eV, respec¬ 
tively. The interaction energy depends also only weakly 


TABLE V. Interaction energies of water molecule with 
NaCl(100) at a high coverage as calculated by different meth¬ 
ods. 


Method 

E int (meV) 

MP2 

-420 

MP2+estimated CC correction 

-430 

EXX+RPA 

-383 

+rSE 

-430 

+GWSE 

-410 

HF+RPA 

-450 


on the k-point sampling, we have used up to 3x3x1 
k-points (the cutoff for the auxiliary basis in this case 
was 250 eV). After accounting for basis set convergence 
and k-point convergence, we obtained an estimate of the 
molecule-surface interaction energy of E- lnt = —420 meV. 
Corrections beyond MP2 are expected to be small. In 
the work of Li et al££, CCSD(T) calculations lead to 
a 10 meV stronger binding energy than MP2. Further 
corrections will arise from the use of hard PAW poten¬ 
tials, for example, the RPA binding increases by 15 meV 
when small core potentials are used, but since we com¬ 
pare MP2 and RPA with similar setups this is irrelevant 
in the present case. 

Our results are summarized in Table As a refer¬ 
ence, we use our MP2 value corrected with the post-MP2 
correction of Li et al., which yields, in total, a binding 
energy of —430 meV. As expected, RPA underestimates 
this value, by about 50 meV. Adding the singles correc¬ 
tion (rSE) leads, in this case, to a perfect agreement with 
the estimated reference data. This is also in agreement 
with the calculations of Ren et al. performed on the 
S22 test set.— They found a very good agreement be¬ 
tween the reference data and RPA+rSE for those dimers 
in S22 that are bound by hydrogen bonds. The hybrid 
HF+RPA overestimates the reference interaction energy, 
which is also in agreement with the findings of Ren et 
al. In this case, the singles in the GW approximation 
(GWSE) do not quite work as well, but the agreement 
with the reference data is still reasonable. 


E. Covalent and metallic solids 

To evaluate the influence of the singles on the lattice 
constants of covalent and metallic solid state systems, we 
show the equilibrium lattice constants for selected mate¬ 
rials in Tab. ED We also use the present opportunity 
to evaluate whether improved PAW potentials have an 
effect on the equilibrium lattice constants. As shown in 
one of our recent work, 7 ' quasiparticle energies as eval¬ 
uated in the GW approximation can have large errors, 
since the PAW projectors possibly do not span the un¬ 
occupied orbital space accurately. As a remedy to this 
problem, we have suggested to use PAW potentials with 
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TABLE VI. Lattice constants (in A) of selected semiconductors, insulators and metals. The zero-point vibration corrected 
experimental lattice constants are taken from Ref. |84|. The table also summarizes the mean relative errors (MRE) and the 


mean absolute relative errors 

(MARE) in percent. 





RPA 

NC-PAW 

RPA 

std-PAW 

+rSE 

std-PAW 

+GWSE 

std-PAW 

exp 

Na 

4.140 

4.195 

4.188 

4.200 

4.214 

Al 

4.022 

4.034 

4.016 

4.020 

4.018 

Rh 

3.808 

3.807 

3.818 

3.807 

3.794 

Pd 

3.895 

3.893 

3.929 

3.893 

3.876 

Cu 

3.601 

3.600 

3.659 

3.605 

3.595 

Ag 

4.051 

4.074 

4.086 

4.073 

4.062 

C 

3.562 

3.571 

3.577 

3.575 

3.553 

Si 

5.428 

5.438 

5.445 

5.445 

5.421 

Ge 

5.623 

5.634 

5.635 

5.632 

5.644 

LiF 

3.994 

3.993 

3.998 

3.993 

3.972 

LiCl 

5.076 

5.079 

5.071 

5.073 

5.070 

NaF 

4.551 

4.607 

4.615 

4.612 

4.581 

NaCl 

5.547 

5.588 

5.576 

5.580 

5.569 

MgO 

4.200 

4.217 

4.230 

4.226 

4.189 

SiC 

4.353 

4.367 

4.374 

4.373 

4.346 

AIN 

4.366 

4.382 

4.388 

4.388 

4.368 

A1P 

5.460 

5.470 

5.474 

5.473 

5.451 

AlAs 

5.646 

5.666 

5.668 

5.663 

5.649 

GaN 

4.493 

4.508 

4.510 

4.508 

4.509 

GaP 

5.442 

5.446 

5.459 

5.462 

5.439 

GaAs 

5.620 

5.643 

5.641 

5.634 

5.640 

InP 

5.869 

5.871 

5.888 

5.887 

5.858 

InAs 

6.036 

6.062 

6.088 

6.089 

6.049 

InSb 

6.471 

6.463 

6.474 

6.451 

6.473 

MRE 

-0.06 % 

0.24 % 

0.44 % 

0.29 % 


MARE 

0.30 % 

0.31 % 

0.55 % 

0.37 % 



norm-conserving partial waves. These unfortunately in¬ 
crease the computational cost, sometimes, even quite 
significantly. In Tab. ED the first column reports the 
lattice constants evaluated using such norm-conserving 
PAW potentials. In the present calculations, to attain 
the highest possible accuracy, the entire lower lying core 
shell was included in the correlated calculations, except 
for oxygen, carbon, nitrogen, fluorine (2 p elements) as 
well as phosphorus and chlorine. For instance for Si and 
Al, the 2s and 2 p states were treated as valence states, 
for In and Sb the 4 d, 4s and 4 p states were fully included. 
The calculations were performed using 6x6x6 k-points 
and 8x8x8 k-points for gapped systems and metals, 
respectively. Increasing the k-point set for selected semi¬ 
conductors and insulators from 6x6x6 to 8x8x8 
changed the lattice constants by less than 0.1 % (C, Si, 
SiC, LiCl). For metals, an increase of the k-point set 
to 10 x 10 x 10 k-points changed the lattice constants 
also only by typically 0.1 % (the results for the transi¬ 
tion metals and std-PAW are reported for 10 x 10 x 10 


k-points). This suggests that the lattice constants are k- 
point converged to about 0.1 %. Errors incurred by the 
finite plane-wave basis set are of the same order, so that 
we estimate the accuracy of the present calculations to 
be about 0.2-0.3 % in the lattice constants (or 1 % in the 
volumes). 

For the RPA the mean relative error with respect to 
the zero-point corrected experimental lattice constants 
is just 0.06 % in the present calculations, whereas the 
mean absolute error is about 0.3 %. We note that this 
is within the estimated error bars of our calculations. 
It is therefore futile to seek for any systematic errors: 
the RPA seems to be able to predict lattice constants 
in almost perfect agreement with experiment. Only for 
the Na metal, the lattice constant is obviously signifi¬ 
cantly underestimated (excluding Na from the calcula¬ 
tions, the MRE and MARE drop to 0.01 % and 0.25 % 
respectively). There are very few density functionals 
that yield a similar accuracy. In fact, the present re¬ 
sults slightly surpass those for the HSEsol functional^ 
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FIG. 2. Energy-volume curves of Pd for RPA, RPA+rSE and 
RPA+GWSE. Clearly the rSE increases the volume. Absolute 
energies are shifted to yield zero at the minimum. 


The most commonly used functional, the PBE functional, 
overestimates the lattice constants by about 1.3 %, and 
even the PBEsol functional yield a mean absolute rela¬ 
tive error of 0.46 % for a slightly larger set”. The use of 
not norm conserving PAW potentials increases the lattice 
constants, on average by 0.3 %. Also the mean absolute 
relative error increases slightly from 0.3 % to 0.35 %. For 
most elements, the differences between standard and NC 
PAW potentials are small. However, they can approach 
up to 0.4 % for elements with 3d and 4d semi-core states 
and (AlAs and GaAs) and up to 1 % for ionic compounds 
with strongly polarizable cores (Na). Note that in ionic 
solids, vdW interactions are sizable, since the Na 2s and 
2 p core electrons are hardly screened and interact via 
van der Waals interactions with the neighboring halide 
atoms. 

The origin for the increase in the lattice constants from 
the NC PAW potentials to the standard PAW poten¬ 
tials is that the standard potentials underestimate the 
polarizability of the core and this in turn yields too large 
lattice constants. For most applications, this small er¬ 
ror of the lattice constants should be acceptable, how¬ 
ever. We finally note that the present values for the 
standard potentials are also in good agreement with our 
first publication”, although the potentials have been im¬ 
proved since our previous calculations published in 2008. 
Specifically, the present set of standard potentials pre¬ 
serves the norm of the pseudo-orbitals better (although 
not perfectly, as the NC PAW potentials do), which 
in turn increases the core polarizability and decreases 
the lattice constants compared to the original values in 

Ref. 0. 

Because the compute cost for the NC PAW potentials 
is very high, we have evaluated the singles only for stan¬ 
dard potentials. The important result is that the sin¬ 
gles hardly change the lattice constants, except for some 
transition metals where an increase by 1 to 2 % is ob¬ 
served for Cu and Pd. Fig. [2] shows that singles indeed 
shift significantly the equilibrium volume to larger val¬ 
ues, clearly worsening the agreement with experiment. 
GWSE rectifies this, and yields almost identical values 



FIG. 3. Measure of the “distance” between HF Green’s 
function and RPA Green’s function. See text for details. 


to the RPA for all considered elements. As we will dis¬ 
cuss in the next paragraph, the RPA density matrix for 
metals is seemingly very different from the HF density 
matrix, and most likely close to the DFT density matrix. 
This suggests that HF singles are not adequate for met¬ 
als, whereas, we expect GW singles to be accurate across 
all systems. 

Fig. [3] indicates how close the GW density matrix is 
to the HF density matrix. In order to measure this, one 
has to introduce a metric to sensibly determine the dif¬ 
ference. An obvious choice is the total energy difference 
between the Hartree-Fock energy evaluated using the HF 
density matrix (rSE) and the Hartree-Fock energy eval¬ 
uated using the RPA/GW density matrix 

Ti'[7rpa/ihf — 7hf^hf]- 

To present the values in a concise manner, we divide this 
by 

Ti'[7dft^hf — 7hf^hf], 

where 7dft is the DFT density matrix and finally take 
the square-root. The reason for including the square- 
root is that the functional is variational and therefore 
quadratic around 7hf, since 7hf is the groundstate den¬ 
sity matrix of /ihf- If the value is 0, the RPA density 
matrix coincides with the HF density matrix, whereas 
for 1 it is closer to the DFT density matrix. One clearly 
recognizes that the RPA density matrix is generally quite 
close to the HF density matrix for light elements and in¬ 
sulators. However, for metallic materials Na, Pd, Cu, 
Rh, and Ag, as well as for small gap semiconductors and 
semiconductors with heavy elements (Ge, Ga and In com¬ 
pounds) this is not the case. The variational properties, 
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discussed above, suggest that an evaluation of the mean 
field contribution using HF orbitals will be generally su¬ 
perior to an evaluation using DFT orbitals. In particular, 
for large gap systems, such as rare gas solids, ice, but also 
C, Si and SiC, LiF and MgO, the differences between rSE 
and GWSE mean held contributions are tiny and only of 
the order of 10 % (essentially the square of the distance 
shown in Fig. 0 . For metals, this is, however, clearly not 
the case, and already shown in the previous paragraph, 
erroneously increases the lattice constants compared to 
RPA or RPA+GWSE. 


V. DISCUSSION AND CONCLUSIONS 

The present work is devoted to the performance of the 
random phase approximation for extended systems if sin¬ 
gles contributions are taken into account. The first part 
of the paper focuses on the derivation of the singles within 
the adiabatic-connection fluctuation-dissipation frame¬ 
work. Not unexpectedly, this derivation yields exactly 
the same contributions as the singles originally suggested 
by Ren and coworkers^ because standard Rayleigh- 
Schrodinger perturbation theory and coupling-constant 
integration are identical. The coupling-constant integra¬ 
tion and the formulation used here has, however, the ad¬ 
vantage that it gives a very clear picture of what the sin¬ 
gles describe. They account for the “mean field” energy 
changes from the non-interacting Kohn-Sham reference 
system to the interacting system, where we define mean 
field as those contributions arising from a changes of the 
one-electron density matrix. This makes it very clear 
why cohesive energies are increased when going from the 
DFT mean held description to the HF mean held descrip¬ 
tion: the latter contracts the orbitals and thus reduces 
the Pauli repulsion between the atoms or molecules. 

Renormalized singles can be also derived in the present 
framework and the hnal equation for them is particu¬ 
larly revealing (compare Equ. (1321) 1. The “renormal¬ 
ized” singles describe the energy difference between the 
Hartree-Fock eigenvalues and the diagonal of the HF ma¬ 
trix evaluated using DFT orbitals (compare Equ. (1531) '). 
This is not exactly identical to the renormalized singles 
suggested by Ren et alM although we found, in practice, 
that our simpler equation gives virtually the same results 
as Ren’s renormalized singles and, as a bonus, it is trivial 
to implement and most likely already available in most 
codes. 

Inspired by the simple physically transparent form of 
the singles, we have also suggested an alternative form 
for the singles that relies on the RPA-density matrix in¬ 
stead of the HF-density matrix (compare Equ. (1351) 1. We 
have termed this correction GWSE, singles in the GW 
approximation. Clearly this description should be supe¬ 
rior to the simple HF description, as the RPA density 
matrix should be fairly close to the exact groundstate 
one-particle density matrix. In practice, for large band 
gap systems, such as rare gas solids, ice and adsorption of 


water on NaCl, differences between the rSE and GWSE 
are small. This suggests that the HF density matrix is 
often astonishingly close to the RPA density matrix for 
systems with light atoms and large band gaps (compare 
Fig. 0. In such cases, the rSE can be used instead 
of the GWSE with little loss of accuracy. We believe 
this explains why the rSE approximation was so success¬ 
ful for molecules. For metals and heavier elements, the 
approximation becomes increasingly worse and an erro¬ 
neousness increase of the lattice constants is observed for 
some transition metals for the rSE approximation, which 
is rectified by the GWSE. 

We finally note that we have used this paper to revise 
our lattice constants for main group elements and some 
transition metals using highly accurate norm-conserving 
potentials. Although, the differences to the previously 
published values are usually small and only of the order 
of 0.3-0.5 %, the present values should be used as future 
reference. 
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Appendix A: Fluctuation-dissipation expression 


Here we briefly derive the expression for the correlation 
energy, Equ. © starting from Equ. (0. The integrand 
in Equ. can be also written as 



^2,A(ri,r 2 ) 

In -r 2 | 


d'Vi d 3 r 2 , 


(Al) 


where the two particle pair density is defined as 
« 2 ,A(ri,r 2 ) = (T A |^ t (r 2 )^ t (r 1 )V’(ri)i/'(r 2 )|T A ), (A2) 


and rp and i/j' are the Fermionic annihilation and cre¬ 
ation operators. Using the common Fermionic anti- 
commutator relations one obtains: 

ti 2 ,A(n,r 2 ) = (4' A |'0 t (r 2 )V’(r 2 )^ t (r 1 )-i/(ri)|T A ) 

- <5(ri - r 2 ) ('F A |'0 t (r 2 )V’(ri)|’F A ) . (A3) 

S. v ✓ 

«i, A (ri;r 2 ) 


The response function is defined by Eq. m, with the 
density fluctuation operator in second quantization given 
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by: Using this expressions for both coordinates and simplify¬ 

ing the expression, the response function takes the form: 

6n A (T,ri) = V’ t (ri,r)^(ri,r) - ra A (ri). (A4) 

I 

X A (r,ri,r 2 ) = ~('I , A |T[V’ t (ri, T)t/j(r 1 , r)^ 1 '(r 2 ,0)V’(r 2 ,0)]|^ A ) + n A (ir)n A (r 2 ). (A5) 

We are interested in r —>• 0~, implying the ri operators act first; after time ordering, we obtain (as \X is bosonic, no 
sign change): 

X' A (0”,ri,r 2 ) = -(^ A |^ t (r 2 ,0)-i/'(r 2 ,0)V’ t (ri,0 _ )V'(ri,0 _ )|^ A ) + ?i A (ri)n A (r 2 ) . (A6) 

The first term on the r.h.s. now equals the first term on the r.h.s. of Equ. (IA3I) . By substituting for this term in 
Eq. (IA31 we obtain 

« 2 ,A(ri,r 2 ) = -XA(0 _ ,ri,r 2 )+ n A (r!)n A (r 2 ) - <f(r! - r 2 )n liA (r!; r 2 ). (A7) 


This corresponds to the terms depending on A in Eq. ©• 
The terms not depending on A are obtained by analo¬ 
gously rewriting (H/q| t)|A&o) i n Equ. ©. 
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